PROCEEDINGS 

OF SCIENCE 



A novel method for evaluating correlation functions 
in lattice hadron spectroscopy* 



o 

(N 

> 
O 



in: 

(N 
> 

oo 

o 



X 



Justin Foley! Chik Him Wong 

Dept. of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA 

John Bulava 

NIC, DESY, Platanenallee 6, 15738, Zeuthen, Germany 

K. Jimmy Juge 

Dept. of Physics, University of the Pacific, Stockton, CA 95211, USA 

David Lenkner, Colin Morningstar 

Dept. of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA 

Michael Peardon 

School of Mathematics, Trinity College, Dublin 2, Ireland 

We describe a new approach for evaluating hadronic correlation functions which combines 
Laplacian-Heaviside quark smearing with a stochastic estimator of quark propagators. This 
method utilizes noise dilution in a new way to reduce the variance in correlators. The efficacy 
of the new algorithm is demonstrated on a number of systems, including disconnected diagrams 
and multi-hadron correlators, on a small lattice where comparisons with the results obtained with 
exactly determined quark propagators are possible. On larger lattice volumes, the use of exact 
propagators becomes prohibitively expensive, while the stochastic method is still computation- 
ally feasible. 
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1. Introduction 



A long-term aim of our collaboration is the first-principles determination of stable hadron 
masses, and resonance energies and decay widths. To date, considerable progress has been made 
in extracting finite-volume stationary-state energies in the isovector meson, kaon, and baryon sec- 
tors |2|, ^[|. These initial studies, which were performed at relatively heavy pion masses and 
moderate lattice volumes, involved interpolating operators specifically designed to couple strongly 
to single-hadron states. However, at lighter pion masses and larger spatial volumes, explicit multi- 
hadron operators must be incorporated into the analysis. Until recently, the inclusion of multi- 
hadron interpolators has been problematic because the resulting two-point functions may involve 
quark lines which begin and end on the sink time slice and source operators on all spatial sites of 
a time slice. Hence, a general treatment of multi-hadrons is not amenable to conventional point- 
to-all quark-propagator techniques. In analyses involving single-particle operators, disconnected 
contributions to flavor-singlet meson correlators are similarly problematic, and little progress has 
been made to date in determining the 1 = meson spectrum from Monte Carlo simulations. 

In this article, we describe a novel method [|[] for evaluating hadronic correlation functions, 
which combines quark-field smearing with a stochastic estimator. The new method facilitates the 
precise evaluation of multi-hadron and flavor-singlet meson correlators, at a significantly lower 
computational cost than previous approaches. First results obtained on intermediate lattice volumes 
are presented, and we briefly discuss more recent studies performed on larger lattices. 



2. An alternate quark-field smearing scheme 

Our approach relies on the fact that generally in a spectroscopy calculation, gauge-covariant 
smearing is applied to the quark fields in hadron interpolators in order to reduce the coupling of 
these operators to very high-lying states. In the resulting correlation functions, the quark propagator 
M _1 is always sandwiched between smearing operators: SM~ l S, where y = S\j/ is the smeared 
quark field. We define Laplacian-Heaviside or LapH, quark-field smearing by 

S = 0(A + a 2 ) , S ab (x,y) « 4 4 ,>, £ v® (x) vf (y)* , (2. 1) 

k=\ 

with 

& ab {x,y)vf (y) = -A,-(x 4 )vi , ' ) (x), v«*(*) v U)(*) = 5« (2.2) 

where A is a gauge-covariant Laplacian operator constructed from stout-smeared [[J link variables 
U : 

&ab(x,y) = t {u?(x)8(x+k,y) + uZ ab (y)8(x-k,y)-25(x,y)8 ab }. (2.3) 
k=i L J 

The use of smeared link variables in the quark-field smearing operator has been shown to signifi- 
cantly reduce the variance in hadronic correlation functions [J7|]. The eigenvalues of the Laplacian 
operator, — A,-, are negative(see Fig.[T]), and the Heaviside function in Eq. 2.1 imposes a cutoff ( 
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Figure 1: The plot on the left shows the Monte Carlo estimate of low-lying A,'s computed on a 16 3 spatial 
lattice on dynamical ensembles with two different pion masses using stout-smeared link variables. The pion- 
mass dependence of the Laplacian spectrum is seen to be very mild. The right-hand plot shows Laplacian 
spectra computed on 2+1 flavors simulations performed on two different lattice volumes. All other simu- 
lation parameters match. The dashed horizontal lines are guides for the eye, showing that the density of 
eigenmodes grows linearly with the spatial lattice volume. 

parametrized by a 2 ) on high-momentum modes. The level of quark smearing is increased by low- 
ering the cutoff, thereby excluding more Laplacian eigenmodes. To evaluate hadronic correlation 
functions involving LapH-smeared quark fields, one does not need to compute the quark-propagator 
components directly, instead only the matrix elements of the propagator between a subset of Lapla- 
cian eigenvectors are required. For a high enough level of quark smearing, i.e., a small enough 
subset of eigenvectors, it is feasible to compute hadron correlation functions which are currently 
beyond the reach of conventional lattice techniques [|5j|. In practice, however, the volume depen- 
dence of the Laplacian eigenmode distribution is a significant limitation when computing SM~ l S 
exactly. Although the quark-mass dependence of the Laplacian spectrum is mild, the density of 
eigenmodes scales linearly with the spatial lattice volume. Therefore, while an exact treatment of 
SM~ l S may work well on smaller lattice volumes as it stands, the increased density of eigenmodes 
makes this approach impractical on larger lattices. 

2.1 Stochastic estimation 

To mitigate the volume dependence, the LapH smearing scheme can be combined with a 
stochastic estimator. To proceed, we write 




(2.4) 



E (tjtj ' 1 denotes an expectation value over the outer product of noise vectors whose components 
satisfy 



E{T)ia(t)) = 0, E(7] ia (0 % (?')) = $AA', ( 2 " 5 ) 



where i(J) and a(j3) are eigenmode and spin indices, respectively. Note that the noise vectors have 
neither color nor spatial indices. The number of quark-matrix inversions required to estimate the 
smeared quark line in this way is determined by the number of stochastic vectors employed rather 
than the number of Laplacian eigenvectors used in the smearing. However, in general, this naive 
estimate for the smeared quark line is too noisy to be of practical use, and we must apply a variance 
reduction scheme to obtain an improved stochastic estimate. 

In a number of systems, noise partitioning, or dilution, turns out to be a particularly effec- 
tive variance reduction technique [B]. Dilution is implemented by choosing noise vectors whose 
components have unit norm, i.e., Z n or £7(1) noise, and partitioning the noise vector indices, eigen- 
mode, time and spin, into disjoint sets. To each set of indices, d, we assign a projection operator 
pi d ] which acts on the noise vectors, with components 

p ulp M = 1 if ( ? ' i '' a ) and ('M G d ' 

P% ] a p[t,t') = otherwise. (2.6) 

The naive noise average E (t]tj M is replaced by an expectation value involving diluted noise vec- 
tors E {r\ ^ r\ ^) , where rj M = P^l r\ . In the maximal dilution limit, the exact smearing scheme 
is recovered. We refer to the combination of the LapH smearing scheme with a dilute stochastic 
estimator as the Stochastic LapH method. 

3. Implementation 

In practice, the estimates for the smeared quark propagator are implemented in terms of pseud- 
ofermion fields which we call quark-line ends. For each diluted noise vector T\ ^ , we form a quark- 
line source p ^ , with components 

pJ^HEnJaWfcM. (3.i) 

i 

For each source field, there is a corresponding quark-line sink <j>^ = SM~ l p^ . The smeared quark 
line can then be written 

SM- i S = '£EU [d] p [d]t ). (3.2) 

In addition, the fact that the quark propagator satisfies 

M- l (y;x) = Y 5 [M- l (x;y)] Jf Y 5 (3.3) 

and the hermiticity of the quark smearing operator lead to a second expression for the smeared 
quark line 

SM- l S = J £E(pWj W ), (3.4) 
where p [d] = fcpM and M = y 5 M. 
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On a typical lattice, the complete set of pseudofermion fields needed to implement the quark- 
line estimate is too large to save to disk. Instead, to store the quark-line estimates on a given 
configuration, we save the set of low-lying Laplacian eigenvectors used in the smearing, the diluted 
noise vectors r\ W 1 , and, for each noise vector, a tensor with elements 

U(i,a,x 4 ;d) = J>£ (3.5) 

x a 

Using these components, the pseudofermion source and sink fields can be quickly reconstructed 
when needed. Only the Laplacian eigenvectors require a significant amount of storage, which is 
however independent of the number of quark lines to be estimated and the number of noise vectors 
and the dilution schemes used. 

To illustrate how the quark-line ends can be combined to form a hadron correlator, consider 
the simple meson two-point function 

c nn (t) = (o|f n m xff(t)f alwio) |o), (3.6) 

where Q. m (Q. n ) is an arbitrary combination of gauge-covariant spatial displacement operators and 
spin matrices acting on the quark fields. For compactness, all indices apart from the time label have 



been suppressed. Wick contracting the quark fields, and using both Eq. 3.2 and Eq. 3.4, we obtain 
an estimate for the connected component of this correlator 



(3.7) 



where the superscripts r and s label different noise sources and (...) denotes an average over the 
gauge ensemble. Other estimates for the connected correlator, involving different combinations of 



'barred' and 'unbarred' quark line ends, are also possible. However, since Eq. 3.7 involves the 



quark sources pf r \ pf s ^ on a single time slice, t = 0, only, it has the advantage that both quark-line 
estimates are automatically fully diluted in time. 

Defining the mesonic 'operators' [ £J s] s] =^ rlt O OT # ] and ~& [ ^\ = (p^ ]t ft„pi* ] ) ^ = 

pf^Q-Xp^ allows the connected correlator to be written in a more compact form 

c:m«-(^L I <A (o<^ (o)). 0.8) 

The factorization of the correlator estimates into contributions from different hadron operators is an 
extremely useful feature of the method. It facilitates the construction of correlator matrices involv- 
ing large sets of interpolating operators and simplifies the evaluation of multi-hadron correlators. 



4. First tests 

Tests of the Stochastic LapH algorithm on connected hadron correlators, such as isovector- 
meson and baryon two-point functions, have shown that, with a judiciously chosen dilution scheme, 

'in practice, we simply record a random initial seed, the dilution scheme, and a rule for generating the noise vector 
from the seed. 
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this method yields statistical errors that are comparable to the errors obtained using an exact treat- 
ment of SM~ X S , at a considerably lower computational cost. Moreover, the new method is more 
efficient than using diluted stochastic estimates with noises introduced on the entire lattice, as de- 
scribed in Ref. This is clearly demonstrated in Fig. ||, which shows the ratios of standard 
deviation on a single time slice of a nucleon correlator evaluated using different stochastic methods 
to that obtained from exact LapH. The solid symbols denote the results obtained by introducing 
noise on the lattice. The open symbols are results from the Stochastic LapH method. The nucleon 
operator in question has gauge-covariant displacements in three directions. These measurements 
were performed on 100 2 + 1 flavor gauge-field configurations with a pion mass of approximately 
400 MeV, on a 20 3 x 128 anisotropic lattice. The spatial lattice spacing was approximately 0. 12 fm, 
and the ratio of spatial to temporal lattice spacings was 3.5. The smearing uses the 64 lowest-lying 
Laplacian eigenmodes on each time slice. The statistical uncertainties on the correlators were es- 
timated using the Jackknife algorithm. The data here correspond to a temporal separation between 
the source and sink operators of t = 5a t , although the same qualitative behavior is observed for 
other temporal separations. The minimal level of dilution included in this figure is full time di- 
lution. The reason is that for a connected correlator evaluated on a single source time slice, the 
inversion cost is the same regardless the time dilution scheme. Any non-full time dilution schemes 
would only impose larger variance from the noises without any reduction of computational costs. 
The eigenmode dilution schemes considered include partitioning the Laplacian eigenmodes into 
blocks of adjacent modes, as well as choosing subsets consisting of eigenmodes that are separated 
by some number n of intermediate modes, known as interlacing. The Stochastic LapH estimates use 
three independent Z4 noise vectors - the minimum number required to obtain unbiased stochastic 
estimates. It is possible to improve the correlator estimates by averaging over different orderings of 
the noise vectors. However, this average has not been performed on the examples presented here. 

For a given number of quark-matrix inversions, the error on the correlator computed using 
Stochastic LapH is much smaller than the error on the correlator estimate involving stochastic 
quark propagators. Fig. |2| also shows that there is no significant difference between using blocked 
eigenmode dilution projectors or interlaced eigenmode projectors. 

Fig. [| compares the same error ratios for Stochastic LapH estimates of the nucleon correlator 
on two different spatial lattice volumes. Results for the 20 3 spatial volume are plotted together with 
results from a smaller, 16 3 lattice. 32 Laplacian eigenvectors are used to smear the quark field on 
the 16 3 volume, such that eigenvalue cutoffs on the two lattices are approximately the same. All 
other input parameters are kept fixed. Increasing the spatial lattice volume for a fixed eigenvalue 
cutoff and dilution scheme does lead to some increase in the statistical error. However, this increase 
is modest for all dilution schemes considered and the difference in errors between the two volumes 
decreases with higher levels of dilution. 

5. Disconnected and multi-hadron correlators 

While the results presented for the connected hadron correlators are promising, the true test 
of the Stochastic LapH method lies in its application to disconnected diagrams and multi-hadron 
correlators. In this section, we compare results obtained for these systems using both the Stochastic 
and exact LapH schemes. Due to the large number of quark-matrix inversions required to imple- 
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Figure 2: Ratio of standard deviation on time slice 5 of a triply-displaced nucleon correlator evaluated 
using different stochastic estimators to the standard deviation on the exact LapH estimate of the correlator 
at the same time separation. The filled symbols denote measurements obtained using noise introduced on 
the entire lattice and the open symbols are results from the Stochastic LapH approach. For every dilution 
scheme considered, Stochastic LapH significantly outperforms the lattice noise approach [pfl. 
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Figure 3: The statistical errors on a triply-displaced nucleon correlator at a source-sink separation of 5a t 
on two different lattice volumes. The correlator was computed using the stochastic smearing scheme with 
varying levels of dilution. The open symbols are results obtained on a 16 3 spatial lattice volume, where the 
quark-smearing operator used 32 Laplacian eigenvectors. Closed symbols denote results on a 20 3 spatial 
volume using 64 Laplacian eigenmodes. The growth in statistical error which comes from increasing the 
spatial volume while the smearing cutoff and the number of quark-matrix inversions per configuration are 
kept fixed is modest [ph. 
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Figure 4: Disconnected contribution to an / = pseudoscalar correlator evaluated on 52 2 + 1 flavor con- 
figurations on a 16 3 x 128 lattice. The pion mass is approximately 400 MeV. The circles are results obtained 
using the exact LapH smearing scheme, and the triangular data points are Stochastic LapH estimates. The 
dilution scheme is interlaced 16 in time, full spin dilution, and interlace 8 in Laplacian-eigenmode space. 
Although the error estimates from both methods are similar, for this particular correlator, the stochastic 
approach requires a factor of 32 fewer quark matrix inversions than are needed in the exact scheme. 

ment the exact scheme, the comparison could only be performed on the smaller 16 3 x 128 lattice on 
a limited number of configurations. The exact method requires additional quark-matrix inversions 
for each sink time slice. However, the number of quark-matrix inversions needed for the stochastic 
estimate is drastically reduced by applying only partial time dilution to the quark lines at the sink. 
In that case, the optimal dilution schemes are expected to involve projectors which are interlaced 
in time. 

Disregarding contact terms, the disconnected contributions to flavor-singlet meson correlators 
can be evaluated using a single diluted noise vector. We have computed the disconnected contribu- 
tions to the isosinglet pseudoscalar and scalar correlators using the interpolating operators ij/iYsYl 
and xffiYi, respectively, where the subscript / denotes a light quark flavor. Fig. |]and Fig. |5] show 
results for the disconnected contributions. These results were obtained on 52 configurations and 
the correlators were averaged over 128 source time slices. The legends list the dilution schemes 
used in the quark-line estimates. [F,F,F] denotes full dilution in time, spin and eigenmode space, 
which corresponds to the exact smearing scheme. [I16,F,I8] indicates that dilution projectors are 
interlaced in time, with projectors containing every sixteenth time slice; the Stochastic LapH es- 
timates employ full spin dilution, which was found to significantly reduce the variance in certain 
disconnected diagrams, and interlace eight eigenmode dilution. In both channels, the exact smear- 
ing scheme result and the dilute Stochastic LapH estimate have similar errors. However, the total 
number of quark-matrix inversions required for the exact estimate is 128 x 4 x 32 = 16384, while 
the stochastic method involves just 16x4x8 = 512 matrix inversions. 

As an aside, the reader may have noticed a discrepancy between the exact and Stochastic LapH 
results for the contact term in Fig. |], and, more obviously, in Fig. ||. This is due to the use of a 
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Figure 5: Disconnected contribution to an isosinglet scalar correlation function estimated using the exact 
and Stochastic LapH methods. The correlator shown was obtained after the subtraction of a large vev term. 
The discrepancy at T = arises because in the Stochastic LapH estimate a single diluted noise vector has 
been used for both the source and sink operators. 

single noise vector for both source and sink operators, giving a biased correlator estimate when the 
temporal separation z = 0, but having no effect on the measurement of spectral quantities. 

Note that the scalar correlator shown in Fig. || was obtained after the subtraction of a large vac- 
uum expectation value contribution. The quality of this signal, obtained on just 52 configurations, 
suggests that considerable progress can be made in the isosinglet meson sector using the Stochastic 
LapH method. 

In the evaluation of disconnected diagrams, it is also possible to use different dilution schemes 
for the source and sink operators. One can, for example, use fully time-diluted noise vectors, 
borrowed from the calculation of the connected correlator component, to estimate the contracted 
source operator in the disconnected contribution. However, in this case, an average over all time 
slices is impractical. We tested this alternate estimate of the disconnected contribution, averaging 
the correlator over four randomly chosen source times. However, on a moderate number of config- 
urations using a lower level of time dilution but averaging correlator estimates over all lattice time 
slices was found to give a cleaner signal, although the discrepancy between the estimates decreases 
on larger ensembles. 

The inclusion of multi-hadron operators is a particular challenge for lattice hadron spec- 
troscopy. Not only are the quark-line diagrams much more complicated, but we also need to 
evaluate correlators involving states whose constitutent hadrons have non-zero momenta. 

The box diagram shown in Fig. ^ contributes to two-pion correlators in the 7 = 0,1 sectors. The 
contribution of this diagram to an S-wave two-pion correlator measured on the 16 3 x 128 lattice is 
displayed in Fig. [7|. The interpolating operator used is simply the product of two zero-momentum 
single-pion interpolators. Once again, the circles denote results obtained using the exact smearing 
scheme, and the triangular data points are results from Stochastic LapH. In the stochastic scheme, 
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Figure 6: The two-pion box diagram and the rho-two-pion diagram discussed in the text 
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Figure 7: The box-diagram contribution to the S-wave, two-pion correlator evaluated with exact and 
Stochastic LapH. These results were again obtained on the 16 3 x 128 lattice. The correlator has been aver- 
aged over 4 well-separated source time slices. The stochastic estimate uses just one diluted noise vector per 
quark line. In the stochastic scheme, quark lines connecting source and sink time slices are fully time dilute, 
while quark lines which begin and end on a single time slice use the interlaced 8 time dilution scheme. Once 
again, the results from the stochastic method are as accurate as results from the exact smearing scheme. 



the quark lines connecting operator source and sink time slices are fully time diluted, while the 
quark lines which begin and end on a single time slice employ the interlaced 16 time dilution 
scheme. The stochastic estimate uses just one diluted noise vector per quark line. In both cases, 
the correlator estimates have been averaged over four source times. 

In truth, Fig. |7| is not a fair comparison of the exact and Stochastic LapH schemes, since if 
it is feasible to apply the exact method, one can place source operators on all lattice time slices. 
However, on larger lattices the exact smearing scheme is prohibitively expensive, and Stochastic 
LapH is the only practical means of evaluating multi-hadron correlation functions. This figure does 
demonstrate the potential for obtaining accurate estimates for multi-hadron correlation functions 
using Stochastic LapH. A comparison of the error bars on the correlator estimates indicates that the 
variance in the Stochastic LapH estimate is dominated by gauge-field fluctuations, and no increase 
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Figure 8: Two-point function correlating a yYkV source operator and a P-wave two-pion interpolator at the 
sink. In this case, the signal from Stochastic LapH is somewhat noisier than the exact method, although still 
reasonable. The exact method becomes computationally intractable on larger lattice volumes. 



in the number of noise vectors or level of dilution is required. 

As a final test of the Stochastic LapH method, we consider the two-point function correlating a 
p -meson operator at the source with and two-pion state at the sink. The 1 = 1 two-pion system is in 
a P-wave state, which requires that the pions have non-zero back-to-back momenta. For the source 
operator, we use the simple quark bilinear: Y7k¥- This operator is correlated with an interpolating 
operator for two pions with one unit of momentum each. As in the previous examples, the stochastic 
estimate uses a single noise vector per quark line. In this case, the Stochastic LapH algorithm gives 
a noticeably noisier result than exact smearing. However, the increase in error is moderate, and the 
stochastic result could be improved at negligible extra cost by averaging the correlator over the two 
possible ways of assigning noise vectors to the quark lines connecting source and sink time slices. 



6. Recent developments 

Motivated by the results from the 16 3 lattice, we have recently begun a systematic study of 
the light hadron spectrum on a 24 3 x 128 lattice 2 , using, in addition to single-hadron operators, 
interpolating operators for two-meson and meson-baryon states. To construct the required two- 
hadron correlators, we need, for each valence quark mass, independent estimates for two quark 
lines that begin and end on a single time slice together with estimates for five quark lines connecting 
source and sink times. We use a single (diluted) noise vector per quark-line estimate, and employ 
the same dilution schemes that were applied to the disconnected and two-pion correlators on the 16 3 
lattice. Connected correlators have been averaged over four randomized source times. The number 
of quark matrix inversions needed to perform measurements using Stochastic LapH smearing is 

2 This corresponds to a spatial volume of approximately (2.9 fm) 3 . 
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then a factor of 34 times smaller than the number of inversions needed to implement the exact 
scheme. 

The results obtained to date are promising. In particular, a comparison of isoscalar meson and 
two-pion correlators measured on the 24 3 spatial volume with results from the 16 3 lattice confirms 
that the mild volume dependence of the Stochastic LapH method in these systems. A detailed 
study of the Stochastic LapH algorithm on the 24 3 x 128 lattice will be presented in a forthcoming 
publication [JTo|] - 



7. Summary and conclusions 

We presented a new algorithm for estimating hadron correlation functions involving smeared 
quark fields. Using Laplacian-Heaviside smearing, the number of quark matrix inversions needed 
to accurately estimate certain hadron correlators can be drastically reduced. However, the strong 
volume dependence of the smearing scheme means that an exact treatment of SM~ l S is impractical 
on larger lattice volumes. We have outlined a way of implementing LapH smearing stochastically, 
using noise dilution to control the variance. Tests involving nucleon correlators show that the 
Stochastic LapH approach is more efficient than using dilute stochastic estimates for the quark 
propagators. Crucially, we have found that Stochastic LapH exhibits a mild volume dependence in 
this sector. The efficacy of this approach for disconnected correlators and multi-hadron correlators 
has also been demonstrated. On a small lattice volume, it was possible to compare the stochastic 
implementation of the LapH smearing scheme with the exact method. We found that Stochastic 
LapH gives results that are close to, or, in many cases, as accurate as the results obtained from the 
exact method, at a fraction of the computational cost. Finally, we noted that recent studies using 
Stochastic LapH performed on a larger lattice confirm excellent volume scaling for disconnected 
and multi-hadron correlators. 

This algorithm is an exciting development for our spectroscopy program. It is now possible 
to include multi-hadron operators in our Monte Carlo measurements, and we can begin to extract 



resonance energies and widths from the measured finite- volume spectra QllL [12J , |13Q . Our results 
indicate that the isosinglet meson spectrum is also now accessible, although in that sector one 
also has to take potential mixing with glueball states into account. The challenges that remain to 
properly dilineate the low-lying hadron spectrum should not be underestimated, but, with continued 
technical and theoretical advances, considerable progress can be made toward this goal in the near 
future. 
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